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Abstract: The paper starts out with a computational technique that allows to 
compare all linear methods for PDE solving that use the same input data. This 
is done by writing them as linear recovery formulas for solution values as linear 
combinations of the input data. Calculating the norm of these reproduction 
formulas on a fixed Sobolev space will then serve as a quality criterion that 
allows a fair comparison of all linear methods with the same inputs, including 
finite-element, finite-difference and meshless local Petrov-Galerkin techniques. 
A number of illustrative examples will be provided. As a byproduct, it turns 
out that a unique error-optimal method exists. It necessarily outperforms any 
other competing technique using the same data, e.g. those just mentioned, and 
it is necessarily meshless, if solutions are written "entirely in terms of nodes" 
(Belytschko et. al. 1996 [4]). On closer inspection, it turns out that it coincides 
with symmetric meshless collocation carried out with the kernel of the Hilbert 
space used for error evaluation, e.g. with the kernel of the Sobolev space used. 
This technique is around since at least 1998, but its optimality properties went 
unnoticed, so far. 

1 Introduction 

For simplicity, consider a standard Dirichlet problem 

Lu = / in O , , 

Bu = g in r := dn ^ ' 

with a linear differential operator L and a linear boundary operator B posed 
on a domain VL. Assume that the amount of available information is fixed and 
limited, namely to values 

f{xj) in points Xj g fi, 1 < .? < ™, fr,\ 

g{yk) in points j/fc G F dil, 1 < k < n. 



Under all methods for solving such problems, we ask for the one with smallest 
error that uses this information and not more. To make this more precise. 
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we fix a single point x € ft and assume existence of a true solution u* of 
the problem. Any method using the above information will produce some 
numerical solution u/^ , and we want to single out methods that make the error 
\u*{x) — um(x')\ small for all problems posed that way. We achieve this by 
looking at error bounds 



where the problems and their solutions are allowed to vary in such a way that the 
true solutions lie in some fixed reproducing kernel Hilbert space iJ of functions 
on ri, e.g. a fixed Sobolev space. In this sense, the constant Cm describes the 
worst-case error behavior of the method M. on all problems with solutions in 
H using the same data. We shall show how to evaluate Cm numerically, and 
this will allow fair comparisons between methods using the same data. A few 
examples including important methods like finite elements, gerneralized finite 
differences, and meshless collocation in various forms are provided at the end. 

But one can also ask for a method M that makes Cm minimal over all linear 
methods using the same data. This problem will be written as one of optimal re- 
covery of functions, and it is proven that the optimal solution exists uniquely. It 
can be calculated explicitly, takes the form of a recovery method in reproducing 
kernel Hilbert spaces, and is a meshless method, if users ask for values of the so- 
lution at nodes. Since it is optimal in the above sense, it outperforms errorwise 
any other competing technique using the same data, may it use finite elements, 
finite differences, or local Petrov-Galerkin techniques. Upon closer inspection, 
it turns out to be nothing new, since it is a special case of symmetric meshless 
collocation, computed using the kernel that was used for error assessment . This 
method relies on Hermite-Birkhoff interpolation in reproducing kernel Hilbert 
spaces [3D] and its application to PDE solving was first analyzed in [TJ [H] . 

2 Recovery Problems 

To keep the presentation simple, we stay with the Dirichlet problem example. 
Any linear method that provides an approximate solution value u{x) for a fixed 
point X and uses exclusively the data 1^ , must necessarily satisfy a formula of 
the form 



whatever the weights rj{x) and Sk{x) are. We call this a direct linear recovery 
of the value u{x) from the data on the right-hand side. The existence of (O 
follows from linearity and the restriction to the admitted data, but in general it 
will need some effort to rewrite a classical method in this form. We come back 
to this in the section on specific methods. 



u*{x) — um{x)\ < C'A4|[w*||ff for all u* S H, 




(3) 
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Clearly, the error in ([3]) is 

rn n 

u*{x) - u{x) = u*{x) - ^ rj{x)/^u*{xj) ~ ^ Sk{x)u*{yk), 

where we now have connected the data functions / and g back to the true 
solution u* . The map 

m n 

ex,r,s ■ u* u*{x) -^rj{x)Au*{x.j) -^Sk{x)u*{yk) (4) 
j=i k=i 

is a linear functional in u*. and if it is bounded on some Hilbert space H, we 
have an error bound 

\u*{x) -u{x)\ < \\ex.r,s\\H'\\u*\\H for all u* e H. (5) 

If the norm ||e^^r,s||ff* is evaluated, it precisely describes the worst-case error 
behavior for all problems with solutions in H, because the inequality is sharp 
by definition of the norm of the functional. In other words, one can express the 
error explicitly in percent of and this is what we shall do in the rest of 

the paper. Furthermore, one can ask for a choice r* and s* of the coefficients 
that make the norm minimal, and this will later lead to an optimal method for 
the given reconstruction. 

Note that we do not consider the way numerical solutions of the considered meth- 
ods are actually calculated. The amount of numerical integration needed for 
weak methods, the error committed in integration subroutines, and any multi- 
grid solver techniques are completely irrelevant. Different numerical strategies 
within the same class of methods, e.g. finite element solvers with different ways 
of integration or with different element families can be compared directly as well. 
The only requirement is that the method is rewritten in direct linear recovery 
form ([3]), and then the total error is evaluated, containing all method-specific 
internal features of the technique. 

3 Error Evaluation 

Fortunately, it is easy to evaluate such norms in reproducing kernel Hilbert 
spaces, and Sobolev spaces are examples. In such spaces, the inner product and 
the kernel K : 51 — > have the properties 

fix) = if, K{x, ■))h for all feH,xen 
K{x,v) = iKix,-),K{y,-))Hiorallx,yen 

and for all linear and continuous functionals X, fJ. £ H* , 

{X,ii)H' = {VK{x,-),iiyK{y,-))H 
= X-fiyK{x,y), 
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where the upper index at the functionals denotes the variable that is acted 
upon. Details are in the literature on kernels, with [19] being a fairly complete 
reference and |15l 116] being open access compilations for teaching purposes. 

Global Sobolev spaces 14^™ (R'') in d dimensions and for order m > d/2 in the 
standard Fourier transform definition have the radial Matern reproducing kernel 

K{x,y) = -—\\x~y\\r'^'K„,_,/,{\\x~yh), x^yeM.' (6) 
1 (m) ' 

with the modified Bessel function of the second kind. Local Sobolev 

spaces W™{^1) for domains C M'' are norm-equivalent to the global spaces, 
as long as domains are non-pathological, i.e. they satisfy a Whitney extension 
property or a have a piecewise smooth boundary with a uniform interior cone 
condition. We shall use the above kernels for evaluation of errors in Sobolev 
space. 

To avoid double sums and to pave the way for generalizations, we shall calculate 
the error norms after rewriting ^ in terms of functionals as 

N 
k=l 

with 

Xj{u) = Lu{xj), 1 < j < m, and Xm+i{u) = u{yi), 1 < i < n, N = m + n (7) 

in our special case. It should be clear that other differential operators and other 
boundary conditions will just change the functionals here. 

Then we get the quadratic form 

N 

Wex.cWH' = K{x,x) ~2^Ci{x)X^K{x,z) 

(8) 

+ Y, c^{x)cj{x)XfX'jK{y,z) 

which can be explicitly evaluated if the kernel is known and the functionals 
are continuous, in particular on all Sobolev spaces on which the functionals are 
continuous. For instance, if L is a second-order elliptic operator, this poses 
the restriction to — 2 > d/2 on the Sobolev order m we can use. But this is 
normal, since we focus on methods that use / = Lu* pointwise, forcing / to be 
in with k > d/2 and thus u* e 1^2™ with to > 2 + d/2. All methods that 
use these data are implicitly making this smoothness assumption, even if their 
users think that they are working in in less regular spaces. Plenty of papers 
and books miss this point. In particular, the standard FEM method is usually 
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formulated and analyzed in low-regularity spaces like or H^, but when it 
comes to implementations using data like ([2]), it needs iJ™ with m > 3 in K^. 

All methods based on the data Ai(w*), . . . , Xk{u*) and brought into the recovery 
form 

N 

u{x) = Ci{x)Xi{u*) 
i=l 

generalizing ^ and (jl} can now be plugged into (|8]) to show how good the 
reproduction quality at x is, since (|5j generalizes accordingly. Note that there 
is no linear system to be solved once the recovery formula is known. But since 
the formula is inserted into the positive definite quadratic form ([5]), a small 
final value will necessarily contain quite some amount of cancellation. The 
computational complexity for norm evaluation is 0(N'^). 

This allows a fair comparison of all such methods using the same data, and we 
shall provide examples below. The comparison can be made pointwise, as we 
saw, but for small problems one can plot the function x ||ea;,c||H* to see where 
a method works badly and needs more data. This gives direct information for 
refinement of the discretization. 

4 Optimal Methods 

Since §1 is a quadratic form with a positive semidefinite matrix, it can be 
minimized. The optimal solution coefficients c*{x) solve the system 

N 

^ cl{x)XyXlK{y, z) = \]K{x, z), 1 < j < iV, (9) 

fc=i 

and the system can be proven to be solvable (20]. Clearly, this choice of coef- 
ficients will then outperform all other competitors error-wise. We shall show 
examples later. The minimal value of ([5]) then is 

N 

\\e.,cA\l' = K{x,x)-Y,c,{xrXlK{x,z)>0 

i=l 

and does not contain any matrix. 

The system ([SJ reveals the nature of this method. Indeed, the functions are 
necessarily linear combinations of the functions XlK{-,z). and application of Af 
to the above system shows the Lagrange property Ai(cfc) = Sik, 1 < i, k < N . 
This means that the Ck are the Lagrange basis for general Hermite-Birkhoff 
interpolation of the given data by the functions Xf K {■ , z)„a,nd this is the well- 
known method of symmetric meshless collocation based on |20| and analyzed 
thoroughly in [51 [7] • The optimality of the technique in the sense of this paper 
should have been known since at least since 1997, but it went unnoticed because 
|15| p. 82, (4.2.2)] was not applied to PDE solving at that time. 
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5 Special Methods 

We now give some details on how to evaluate recovery errors for special PDE 
solution techniques. This will be useful for the examples in the final section, 
and we specialize to L ~ —A here. 

5.1 Finite Elements 

The simplest possible 2D finite-element code uses piecewise linear elements on 
the triangles of a triangulation of a domain with piecewise linear boundary, 
and requires /-values only at the barycenters of the triangles. These are the 
points Xj in ^ in the FEM version that we denote by FEMBary below. But 
since usually there are more triangles than nodes, one can also prescribe /- 
values at all interior and boundary nodes to calculate approximate values at 
the barycenters. This usually needs less /-values at the same order of accuracy. 
We call this method FEMNode below. These two FEM variations have different 
recovery formulas ^ and different error functionals (jj} to be compared. 

To get the recovery formulas in the form needed here, users will have to check 
carefully what their FEM code does. We used the MATLAB pdetool setting, 
which does the following. It uses all N triangle vertices zi,...,zjv for set- 
ting up the standard piecewiese linear test and trial functions vi. . . . ,vn with 
Vj{zk) — 5jk and builds the N x N stiffness matrix with entries (Vu^, Vuj)L2 
in the usual way. The right-hand sides (/, i'j)l2i ^ ^ j ^ N are certain linear 
combinations of /-values either at triangle barycenters or at nodes of the neigh- 
boring triangles. These linear combinations form a sparse integration matrix B 
with entries bjk such that the stiffness system without boundary conditions but 
with some form of numerical integration would be 

N N 

Y,i^V^.^VJ)LMz^) = bjkfixk), 1<J<N. 

i=l k=l 

The FEMBary and FEMNode variations have different B matrices and use / 
at different points Xk, but the stiffness matrices are the same. 

The unknowns u{zi), i E D C {1, . . . , N} are known Dirichlet values, and thus 

N 

i^D fc=l ieD 

is the system to be actually solved. Since the u{zi) with i E D are g- values, it 
has the matrix-vector form 

Au = Bf + Cg (10) 
under adequate notation, and each row of 

u = A-^Bf + A^Cg 
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provides one instance of Q for each of the points x = Xi, i ^ D. If we focus on 
the origin as one of these points, we need just one row. 

If users want the error at a non-nodal point x, they have to add a piece wise 
linear interpolation on a triangle containing x, and then the discrete reconstruc- 
tion formula is a linear combination of three rows of the above system. This 
illustrates how numerical integration and interpolatory post-processing both 
enter explicitly into our version of a complete and explicit FEM error analysis. 

5.2 Symmetric Kernel— Based Collocation 

These work on the two point sets X and Y from ([2]) and use linear combinations 

?n n 

u{x) =Y,CjAK{x,x^) +J2dkK{x,yk) (11) 

3 = 1 k=l 

of basis functions derived from a smooth kernel K . The argument in Section 
m shows that these methods yield optimal errors in the Hilbert spaces in which 
their kernels are reproducing. In case of Sobolev spaces, we thus get the error- 
optimal methods this way. 

The numerical process collocates these trial functions at Dirichlet and PDE 
nodes, forming the block system 

m n 

Y^CjA-AyK{x,,x,)+Y,dk^K{x,,yk) = f{x,),l<i<m, 

3 = 1 k=l 
Tn n 

^CjAK{y^,Xj) +^dkK{y,,yk) 9{yi),l<i<n. 
3 = 1 k=i 

The inverse of the coefficient matrix recovers the coefficients Cj and dk from the 
/ and g data, and the numerical solution at x is just a linear combination pip 
of those coefficients. Consequently, the discrete recovery ^ is furnished by the 
inverse of the above "stiffness" matrix, premultiplied by the row vector of the 
kernel values in pip . 

Of course, one can use a special kernel K to calculate the discrete recovery, 
and then use another kernel, e.g. one generating a Sobolev space, for error 
evaluation on that space. We shall do this in the final section. 

5.3 Unsymmetric Kernel— Based Collocation 

In contrast to the previous section, this class of methods takes an additional set 
Z = {zi, . . . , zat} of usually N = m +n points and works on linear combinations 



N 

u{x) = ^ dkK{x, Zk) 
fe=i 



(12) 
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of basis functions, while still using collocation in the point sets X and Y from 
([2]) • This method dates back to early papers of Ed Kansa |10l [TT] and was called 
MLSQ2 as a variation of the Meshless Local Petrov-Galerkin method |2] [T] of 
S.N. Atluri and collaborators. The linear system for the coefficients now is 

n 

^dkAK{xi,Zk) = f{xi),l<i<m, 
fe=i 

n 
k=l 

and the inverse of the coefficient matrix (if it exists, see [9]), premultiplied by 
the vector of values K{x,Zk) of (|12p will yield a row vector for the discrete 
recovery formula ([3]) . If is chosen larger than m + n to increase stability, a 
pseudoinverse of the system coefficient matrix can replace the inverse. This was 
done in the examples of the final section. 

5.4 Meshless Lagrange Methods 

Here, a set Z — {zi, . . . ^ zn} of trial nodes is chosen, and there are shape 
functions ui, . . . , ujv such that trial functions 

N 

"(^) y^^ukix)uizk) 

k=l 

can be written "entirely in terms of nodes". Usually, this implies Lagrange 
conditions Uj{zk) = Sjk, and in many cases the shape functions are defined 
via Moving Least Squares. We do not care here for details, and allow such 
techniques to come in weak or strong form. The strong case collocates in sets 
X and Y like above, forming a system 

JV 

^Auk{xj)u{zk) = f{xj), l<j<m, 

k=l 

JV 

^Uk{yt)u{zk) = g{yi),l<i<n 

k=l 

and the weights of the discrete recovery at Zk will be a row of the preudoinverse 
of the coefficient matrix of this system. 

The weak cases form stiffness matrices and right-hand sides like in the FEM 
situation, and then we get the coefficients of the discrete recovery in the same 
way, involving a special matrix B caring for the numerical integration. 

5.5 Generahzed Finite— Difference Methods 

Here, there are no trial functions, but everything is still expressed in terms of 
values at nodes Z = {zi, . . . , z^v}. In the strong situation, PDE operator values 
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are approximated by formulas like 



N 



(13) 



with localized weights ajk, and this can be done with minimal error in a re- 
producing kernel Hilbert space using the logic of section ID See |6j [17] for more 
details. The linear system then is 



if we use a subset D of Dirichlet nodes like in the FEM case. This is of the form 
pop and we already know how to derive the recovery formulas in such a case. It is 
interesting to see that the matrix with coefficients ajk plays the role of a stiffness 
matrix here, and it is the place where sparsity can be implemented to yield 
local methods with sparse matrices. This was very successfully done in various 
apphcation papers of C.S. Chen, B. Sarler, and G.M. Yao [111 US HH [23l [18] . 
We shall present a simple numerical example in the final section. 

The Direct Meshless Local Petrov Galerkin methods of [121 [S] are weak cases 
of this approach, with other functionals than (|13p being directly approximated 
in terms of values at nodes, and with integrations involving B matrices again, 
like in all other weak methods. We leave this for a future paper dealing with all 
variants of Atluri's Meshless Local Petrov Galerkin technique, and comparing 
them to FEM and kernel methods, optimal or not, sparse or not. 

6 Numerical Examples 

To avoid overloading the paper, we present a simple series of examples. They all 
work on the unit circle for simplicity, and in order to include finite elements, we 
have to use triangulations, even if they are not needed for meshless methods. We 
start with the standard discretization of the unit circle into 8 triangles meeting 
at the origin, followed by three standard finite-element refinement steps halving 
the edges. The problem ([T]) is posed with L = —A, and Dirichlet boundary 
values are always provided in the boundary nodes, which are the yk in In 
all cases, the non-boundary vertices of the triangles are the nodes where we 
want to know the solution, but in the sense of (H} and for simplicity we only 
evaluate the recovery error at the origin x = 0. 

As the FEM variations show, one can work with / values either in barycenters 
of triangles or in vertices of the triangulation. We shall evaluate all examples 
in both situations, denoting the methods by either *Bary or *Node. Details on 
the discretizations are in Table [1] 



N 
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n: number of Dirichlet boundary data points for g values, 

fTiBary'- number of triangles and barycentric data points for / values, 

iTT-Node'- number of vertices and vertex data points for / values, including 
the n Dirichlet boundary vertices, 

DOF: degrees of freedom = number of unknowns ~ rriNode ~ n, 

h: fill distance in the sense of kernel discretizations, describing the maxi- 
mal distance of an arbitrary point of the domain to one of the vertices. 



Case 


n 




ITlNode 


DOF 


h 


CO 


8 


8 


9 


1 


0.2706 


CI 


16 


32 


25 


9 


0.1515 


C2 


32 


128 


81 


49 


0.0768 


C3 


64 


512 


289 


225 


0.0389 


C4 


128 


2048 


1089 


961 


0.0197 



Table 1: Discretization data for the examples on the unit disk 

The plots in Figures [1] to [5] show the errors in Sobolev space of order 3 to 7, with 
order 3 being somewhat out of bounds, and the data for the fine discretization C4 
being polluted by ill-conditioning in various cases. In each figure, the Sobolev 
space for error evaluation is fixed, but the methods might use other kernels. 
The methods are 

FEMBary: FEM with / data in barycenters 

FEMNode: FEM with / data in nodes 

KansaBary: Unsymmetric collocation with / data in barycenters, using 
the order 7 Sobolev kernel 

KansaNode: same with node data 

HOBary: symmetric high-order collocation with / data in barycenters, 
using the order 7 Sobolev kernel 

HONode: same with data in nodes. These two coincide with the optimal 
methods, if evaluated on Sobolev space of order 7, see Figure [S] 

OptBary: optimal method in Sobolev space used for error evaluation, with 
/ data in barycenters 

OptNode: same with data in nodes 

LocNode: a bandwidth 15 method like in section 15.51 using the order 7 
Sobolev kernel, data in nodes. 
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There are serious instabilities that need explanation. They could have been 
avoided in all cases by choosing a smaller kernel scale, but this would have 
changed the Sobolev space norm and concealed the instabilities. Consequently, 
the kernel scale for Sobolev space error evaluation was fixed at 1.0 through- 
out. This does not seriously afi^ect the error evaluation, because the latter just 
consists of a calculation of a quadratic form. 

But it affects the calculation of recovery formulas by kernel methods, including 
the optimal ones. In particular, increasing the smoothness of the construction 
kernel will increase the instability, if no precautions like preconditioning [3l |5] 
are taken. We chose the Kansa* and the HO* kernel methods to work with the 
kernel of Sobolev space of order 7 at scale 1 , but smaller scales would have been 
more stable. Users can easily try different kernels and scales for construction and 
then evaluate in a fixed Sobolev space at a fixed scale to see which construction 
scale gives best results. Example is missing 

The numerical results support quantitatively what is known from experience and 
partly supported by theory. The piecewise linear FEM technique is adapted 
to low-regularity problems, and it shows its second-order convergence in all 
applicable Sobolev spaces from order 4 on. It is clearly inferior error-wise to 
the error-optimal method, which is symmetric collocation in Sobolev space, 
and the difference gets larger for increasing Sobolev order, because the optimal 
method behaves like a p-FEM and increases its convergence rate automatically 
with the Sobolev order. Unfortunately, it suffers from severe ill-conditioning if 
no precautions are taken, and it does not allow sparsity. However, it serves as 
a standard reference to evaluate error performances of all other methods using 
the same data. 

We now list some observations concerning comparisons. The optimal method 
HO* for a fixed high Sobolev order performs well also for lower Sobolev order. It 
adapts automatically to lower regularity. Other competing methods, like Kansa* 
for unsymmetric collocation, behave in the same way. If sparsity is enforced by 
localizing meshless collocation [TH |^ [2TJ [531 HH] , one gets competitive methods 
Loc* that share sparsity with FEM techniques, but also yield high convergence 
orders depending on the bandwidth chosen. 

7 Conclusion and Outlook 

The paper provides a tool that allows an explicit and fully computational as- 
sessment of the error behavior of all linear solvers for all linear PDE problems. 
The exact solution can be unknown, and the error is expressed as a factor of 
the (unknown) Sobolev norm of the true solution. This tool should be applied 
in many more circumstances, e.g. on special and awkward domains, for more 
general differential operators including those of Computational Mechanics, and 
for many other linear solvers, e.g. MLPG methods and boundary-oriented tech- 
niques like the DRM. Any application-oriented paper can, in principle, apply 
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this technique and thus provide a strict worst-case error bound in terms of the 
Sobolev norm of the true solution. Examples of single cases with known solu- 
tions can never be completely satisfactory, but they are the usual practice in 
application papers. 

There is a method that always realizes the optimal error, but it needs further 
work towards numerical stabilization. All other methods should be compared 
to it. and it is an interesting research challenge to see how close one can come 
to the optimal method under sparsity restrictions. 



Errors in Sobolev space of order 3 
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Figure 1: Errors for Sobolev order 3 
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CO 


CI 




C2 




C3 


C4 


h 


2.706e-001 


1.515e-001 


7.682e- 


002 


3.895e- 


002 


1.965e-002 


FEMBary 


6.011e-001 


3.415e-001 


1.737e- 


001 


8.635e- 


002 


4.272e-002 


FEMNode 


6.912e-001 


4.457e-001 


2.411e- 


001 


1.220e- 


001 


6.063e-002 


KansaBary 


3.150e-001 


1.422e+000 


5.888e- 


001 


9.385e- 


001 


6.421e-001 


KansaNode 


3.150e-001 


9.847e-001 


5.131e- 


001 


4.303e- 


001 


6.696e-001 


HOBary 


5.510e-001 


3.538e-001 


1.802e- 


001 


1.245e- 


001 


3.605e+002 


HONode 


1.060e+000 


4.918e-001 


2.498e- 


001 


1.371e- 


001 


3.475e+000 


OptBary 


1.390e-001 


1.309e-001 


1.090e- 


001 


7.280e- 


002 


4.051e-002 


OptNode 


1.391e-001 


1.353e-001 


1.218e- 


001 


9.145e- 


002 


5.508e-002 


LocNode 


5.236e-001 


4.734e-001 


2.591e- 


001 


1.264e- 


001 


6.196e-002 



Table 2: Raw data for Sobolev order 3 







CO 




Cl 




C2 




C3 




C4 


h 


2.706e- 


001 


1.515e- 


001 


7.682e- 


002 


3.895e- 


002 


1.965e- 


002 


FEMBary 


2.223e- 


002 


5.296e- 


003 


1.301e- 


003 


3.494e- 


004 


9.856e- 


005 


FEMNode 


2.547e- 


002 


7.813e- 


003 


1.984e- 


003 


4.844e- 


004 


1.199e- 


004 


KansaBary 


1.995e- 


001 


5.219e- 


002 


7.806e- 


003 


6.309e- 


003 


1.049e- 


002 


KansaNode 


1.995e- 


001 


4.215e- 


002 


6.083e- 


003 


3.478e- 


003 


l.lOle- 


002 


HOBary 


2.173e- 


002 


4.353e- 


003 


9.958e- 


004 


4.171e- 


004 


6.565e- 


001 


HONode 


2.446e- 


002 


5.855e- 


003 


1.479e- 


003 


4.436e- 


004 


9.688e- 


003 


OptBary 


2.163e- 


002 


4.308e- 


003 


9.355e- 


004 


2.174e- 


004 


5.287e- 


005 


OptNode 


2.029e- 


002 


5.770e- 


003 


1.463e- 


003 


3.653e- 


004 


9.117e- 


005 


LocNode 


4.265e- 


002 


7.420e- 


003 


1.651e- 


003 


4.069e- 


004 


1.364e- 


004 



Table 3: Raw data for Sobolev order 4 







CO 




Cl 




C2 




C3 




C4 


h 


2.706e- 


001 


1.515e- 


001 


7.682e- 


002 


3.895e- 


002 


1.965e- 


002 


FEMBary 


e.OlOe- 


003 


1.762e- 


003 


5.658e- 


004 


1.775e- 


004 


5.376e- 


005 


FEMNode 


8.042e- 


003 


2.134e- 


003 


5.451e- 


004 


1.477e- 


004 


4.202e- 


005 


KansaBary 


1.447e- 


001 


1.522e- 


002 


6.089e- 


004 


2.612e- 


004 


1.005e- 


003 


KansaNode 


1.447e- 


001 


1.345e- 


002 


4.163e- 


004 


1.590e- 


004 


1.092e- 


003 


HOBary 


5.223e- 


003 


3.307e- 


004 


3.299e- 


005 


8.055e- 


006 


5.733e- 


003 


HONode 


3.723e- 


003 


4.258e- 


004 


4.881e- 


005 


9.465e- 


006 


1.303e- 


004 


OptBary 


5.031e- 


003 


3.304e- 


004 


3.106e- 


005 


3.298e- 


006 


3.873e- 


007 


OptNode 


3.278e- 


003 


3.977e- 


004 


4.798e- 


005 


5.890e- 


006 


7.463e- 


007 


LocNode 


1.591e- 


002 


8.987e- 


004 


8.243e- 


005 


1.742e- 


005 


8.366e- 


005 



Table 4: Raw data for Sobolev order 5 



REFERENCES 





CO 


CI 


C2 


C3 




C4 


h 


2.706e-001 


1.515e-001 


7.682e-002 


3.895e-002 


1.965e- 


002 


FEMBary 


3.167e-003 


1.219e-003 


4.051e-004 


1.262e-004 


3.780e- 


005 


FEMNode 


3.485e-003 


1.054e-003 


3.226e-004 


9.917e-005 


3.003e- 


005 


KansaBary 


1.135e-001 


6.470e-003 


6.036e-005 


1.451e-005 


1.186e- 


004 


KansaNode 


1.135e-001 


5.501e-003 


3.513e-005 


9.405e-006 


1.349e- 


004 


HOBary 


1.641e-003 


2.955e-005 


1.262e-006 


2.034e-007 


5.496e- 


005 


HONode 


7.159e-004 


4.206e-005 


1.983e-006 


3.158e-007 


1.965e- 


006 


OptBary 


1.600e-003 


2.949e-005 


1.222e-006 


5.964e-008 


1.566e- 


007 


OptNode 


6.730e-004 


3.811e-005 


1.901e-006 


l.llOe-007 


5.458e- 


006 


LocNode 


7.056e-003 


1.521e-004 


1.247e-005 


1.342e-005 


7.800e- 


005 



Table 5: Raw data for Sobolev order 6 





CO 


CI 


C2 


C3 




C4 


h 


2.706e-001 


1.515e-001 


7.682e-002 


3.895e-002 


1.965e- 


002 


FEMBary 


2.337e-003 


9.327e-004 


3.075e-004 


9.490e-005 


2.820e- 


005 


FEMNode 


1.975e-003 


7.244e-004 


2.425e-004 


7.684e-005 


2.339e- 


005 


KansaBary 


9.332e-002 


3.506e-003 


8.576e-006 


1.232e-006 


1.792e- 


005 


KansaNode 


9.332e-002 


2.743e-003 


4.428e-006 


8.672e-007 


2.207e- 


005 


HOBary 


6.633e-004 


2.916e-006 


5.232e-008 


2.147e-009 


2.658e- 


007 


HONode 


1.722e-004 


6.116e-006 


1.202e-007 


1.390e-008 


5.223e- 


008 


OptBary 


6.633e-004 


2.916e-006 


5.232e-008 


2.147e-009 


2.658e- 


007 


OptNode 


1.722e-004 


6.116e-006 


1.202e-007 


1.390e-008 


5.224e- 


008 


LocNode 


3.590e-003 


3.179e-005 


4.802e-006 


1.273e-005 


7.360e- 


005 



Table 6: Raw data for Sobolev order 7 



